Temas

  1. Lectura de archivos de diferentes formatos
    • Base de datos (data frame)
    • Multipolígono
    • Raster
  2. Paquete sf
    • Proyección y reproyección
    • Análisis
  3. Visualización interactiva con mapview
  4. Visualización para publicación con ggplot2

Esta presentación y los datos utilizados está en https://github.com/PPaccioretti/VisualizacionDatosEspaciales

Lectura de archivos

Carga e instalación de librerías para los procedimientos

libs <- c("sf","raster", "ggplot2","mapview")

new.packages <- libs[!(libs %in% rownames(installed.packages()))]
if(length(new.packages)) install.packages(new.packages)

invisible(sapply(libs, library,character.only = T, quietly=T))

Lectura de archivo Raster

DEM<-raster("../Datos/dtm_elevation_merit.dem_m_250m_s0..0cm_2017_v1.0.tif")
DEM
## class       : RasterLayer 
## dimensions  : 2838, 2286, 6487668  (nrow, ncol, ncell)
## resolution  : 0.002083333, 0.002083333  (x, y)
## extent      : -66.12502, -61.36252, -35.19665, -29.28415  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Users\Pablo\Documents\GitHub\VisualizacionDatosEspaciales\Datos\dtm_elevation_merit.dem_m_250m_s0..0cm_2017_v1.0.tif 
## names       : dtm_elevation_merit.dem_m_250m_s0..0cm_2017_v1.0 
## values      : 38, 2749  (min, max)

Lectura de archivo .shp

Departamentos<-st_read("../Datos/deptos_cba_cg")
## Reading layer `deptos_cba_cg' from data source `C:\Users\Pablo\Documents\GitHub\VisualizacionDatosEspaciales\Datos\deptos_cba_cg' using driver `ESRI Shapefile'
## Simple feature collection with 26 features and 14 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -65.77198 ymin: -35.00013 xmax: -61.77089 ymax: -29.50042
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

Lectura de archivo de texto

Suelos<-read.table("../Datos/suelos.txt", header=T)

Cambio a clase espacial (sf). Es importante indicar cual es la proyección de la base de datos mediante el argumento crs =. Si no se le especifica no tomará ningún sistema de coordenadas de referencia.

Suelossf<-st_as_sf(Suelos,coords = c("Xt","Yt"),  crs = 32720)

En que proyecciones están estos objetos?

st_crs(Departamentos)
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(Suelossf)
## Coordinate Reference System:
##   EPSG: 32720 
##   proj4string: "+proj=utm +zone=20 +south +datum=WGS84 +units=m +no_defs"
st_crs(DEM)
## Coordinate Reference System: NA

Reproyectemos el sistema de coordenadas de referencia

st_crs(Departamentos) <- 4326

Para reproyectar el objeto Departamentos podemos tomar la proyección del objeto Suelossf como referencia

Departamentos<-st_transform(Departamentos,st_crs(Suelossf))

Paquete sf

El paquete sf permite un análisis y manejo de archivos espaciales.

Simple features es una manera estandarizada de codificar en computadoras datos vectoriales (puntos, lineas, polígonos). El paquete sf implementa simple features en R y tiene la misma capacidad para datos vectoriales como los paquetes sp, rgeos y rgdal (Pebesma 2018).

Visualización interactiva con mapview

## S4 method for signature 'sf'
mapView(x, map = NULL, pane = "auto",
  canvas = useCanvas(x), viewer.suppress = canvas, zcol = NULL,
  burst = FALSE, color = mapviewGetOption("vector.palette"),
  col.regions = mapviewGetOption("vector.palette"), at = NULL,
  na.color = mapviewGetOption("na.color"), cex = 6,
  lwd = lineWidth(x), alpha = 0.9, alpha.regions = regionOpacity(x),
  na.alpha = regionOpacity(x), map.types = NULL,
  verbose = mapviewGetOption("verbose"), popup = popupTable(x),
  layer.name = NULL, label = makeLabels(x, zcol),
  legend = mapviewGetOption("legend"), legend.opacity = 1,
  homebutton = TRUE, native.crs = FALSE,
  highlight = mapviewHighlightOptions(x, alpha.regions, alpha, lwd),
  maxpoints = getMaxFeatures(x), ...)

Se puede ver una capa:

mapview(Departamentos)
Colorear por valores numéricos de una columna
r mapview(Departamentos, zcol='st_area_sh')

Colorear por valores categóricos de una columna

mapview(Departamentos, zcol='departa')
Diferentes tipos de mapas de base
r mapview(Departamentos, zcol='departa', map.types = c("OpenStreetMap","CartoDB.DarkMatter"))

Seleccionar las categorias que deseamos ver

mapview(Suelossf, legend=TRUE, cex="pH")

Multiples capas en un solo mapa

mapview(DEM) + Suelossf
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 6487668 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  6487668 '

Visualización para publicación con ggplot

ggplot(Suelossf) + geom_sf()

ggplot(Departamentos) + geom_sf(aes(fill=st_area_sh))

Cisualización de objeto de clase RasterLayer

spplot(DEM)

Referencias